// Magic Software, Inc.
// http://www.magic-software.com
// Copyright (c) 2000, All Rights Reserved
//
// Source code from Magic Software is supplied under the terms of a license
// agreement and may not be copied or disclosed except in accordance with the
// terms of that agreement.  The various license agreements may be found at
// the Magic Software web site.  This file is subject to the license
//
// FREE SOURCE CODE
// http://www.magic-software.com/License/free.pdf

#ifndef MGCMATRIX3_H
#define MGCMATRIX3_H

#include "MgcVector3.h"

// NOTE.  The (x,y,z) coordinate system is assumed to be right-handed.
// Coordinate axis rotation matrices are of the form
//   RX =    1       0       0
//           0     cos(t) -sin(t)
//           0     sin(t)  cos(t)
// where t > 0 indicates a counterclockwise rotation in the yz-plane
//   RY =  cos(t)    0     sin(t)
//           0       1       0
//        -sin(t)    0     cos(t)
// where t > 0 indicates a counterclockwise rotation in the zx-plane
//   RZ =  cos(t) -sin(t)    0
//         sin(t)  cos(t)    0
//           0       0       1
// where t > 0 indicates a counterclockwise rotation in the xy-plane.


class MgcMatrix3
{
public:
    // construction
    MgcMatrix3 ();
    MgcMatrix3 (const MgcReal aafEntry[3][3]);
    MgcMatrix3 (const MgcMatrix3& rkMatrix);
    MgcMatrix3 (MgcReal fEntry00, MgcReal fEntry01, MgcReal fEntry02,
                MgcReal fEntry10, MgcReal fEntry11, MgcReal fEntry12,
                MgcReal fEntry20, MgcReal fEntry21, MgcReal fEntry22);

    // member access, allows use of construct mat[r][c]
    MgcReal* operator[] (int iRow) const;
    operator MgcReal* ();
    MgcVector3 GetColumn (int iCol) const;

    // assignment and comparison
    MgcMatrix3& operator= (const MgcMatrix3& rkMatrix);
    bool operator== (const MgcMatrix3& rkMatrix) const;
    bool operator!= (const MgcMatrix3& rkMatrix) const;

    // arithmetic operations
    MgcMatrix3 operator+ (const MgcMatrix3& rkMatrix) const;
    MgcMatrix3 operator- (const MgcMatrix3& rkMatrix) const;
    MgcMatrix3 operator* (const MgcMatrix3& rkMatrix) const;
    MgcMatrix3 operator- () const;

    // matrix * vector [3x3 * 3x1 = 3x1]
    MgcVector3 operator* (const MgcVector3& rkVector) const;

    // vector * matrix [1x3 * 3x3 = 1x3]
    friend MgcVector3 operator* (const MgcVector3& rkVector,
        const MgcMatrix3& rkMatrix);

    // matrix * scalar
    MgcMatrix3 operator* (MgcReal fScalar) const;

    // scalar * matrix
    friend MgcMatrix3 operator* (MgcReal fScalar, const MgcMatrix3& rkMatrix);

    // utilities
    MgcMatrix3 Transpose () const;
    bool Inverse (MgcMatrix3& rkInverse, MgcReal fTolerance = 1e-06) const;
    MgcMatrix3 Inverse (MgcReal fTolerance = 1e-06) const;
    MgcReal Determinant () const;

    // singular value decomposition
    void SingularValueDecomposition (MgcMatrix3& rkL, MgcVector3& rkS,
        MgcMatrix3& rkR) const;
    void SingularValueComposition (const MgcMatrix3& rkL,
        const MgcVector3& rkS, const MgcMatrix3& rkR);

    // Gram-Schmidt orthonormalization (applied to columns of rotation matrix)
    void Orthonormalize ();

    // orthogonal Q, diagonal D, upper triangular U stored as (u01,u02,u12)
    void QDUDecomposition (MgcMatrix3& rkQ, MgcVector3& rkD,
        MgcVector3& rkU) const;

    MgcReal SpectralNorm () const;

    // matrix must be orthonormal
    void ToAxisAngle (MgcVector3& rkAxis, MgcReal& rfRadians) const;
    void FromAxisAngle (const MgcVector3& rkAxis, MgcReal fRadians);

    // The matrix must be orthonormal.  The decomposition is yaw*pitch*roll
    // where yaw is rotation about the Up vector, pitch is rotation about the
    // Right axis, and roll is rotation about the Direction axis.
    bool ToEulerAnglesXYZ (float& rfYAngle, float& rfPAngle,
        float& rfRAngle) const;
    bool ToEulerAnglesXZY (float& rfYAngle, float& rfPAngle,
        float& rfRAngle) const;
    bool ToEulerAnglesYXZ (float& rfYAngle, float& rfPAngle,
        float& rfRAngle) const;
    bool ToEulerAnglesYZX (float& rfYAngle, float& rfPAngle,
        float& rfRAngle) const;
    bool ToEulerAnglesZXY (float& rfYAngle, float& rfPAngle,
        float& rfRAngle) const;
    bool ToEulerAnglesZYX (float& rfYAngle, float& rfPAngle,
        float& rfRAngle) const;
    void FromEulerAnglesXYZ (float fYAngle, float fPAngle, float fRAngle);
    void FromEulerAnglesXZY (float fYAngle, float fPAngle, float fRAngle);
    void FromEulerAnglesYXZ (float fYAngle, float fPAngle, float fRAngle);
    void FromEulerAnglesYZX (float fYAngle, float fPAngle, float fRAngle);
    void FromEulerAnglesZXY (float fYAngle, float fPAngle, float fRAngle);
    void FromEulerAnglesZYX (float fYAngle, float fPAngle, float fRAngle);

    // eigensolver, matrix must be symmetric
    void EigenSolveSymmetric (MgcReal afEigenvalue[3],
        MgcVector3 akEigenvector[3]) const;

    static void TensorProduct (const MgcVector3& rkU, const MgcVector3& rkV,
        MgcMatrix3& rkProduct);

    static const MgcReal EPSILON;
    static const MgcMatrix3 ZERO;
    static const MgcMatrix3 IDENTITY;

protected:
    // support for eigensolver
    void Tridiagonal (MgcReal afDiag[3], MgcReal afSubDiag[3]);
    bool QLAlgorithm (MgcReal afDiag[3], MgcReal afSubDiag[3]);

    // support for singular value decomposition
    static const MgcReal ms_fSvdEpsilon;
    static const int ms_iSvdMaxIterations;
    static void Bidiagonalize (MgcMatrix3& kA, MgcMatrix3& kL,
        MgcMatrix3& kR);
    static void GolubKahanStep (MgcMatrix3& kA, MgcMatrix3& kL,
        MgcMatrix3& kR);

    // support for spectral norm
    static MgcReal MaxCubicRoot (MgcReal afCoeff[3]);

    MgcReal m_aafEntry[3][3];
};

#endif
